Surface wave tomography using dense 3D data around the Scrovegni Chapel in Padua, Italy

A dense single-node 3D seismic survey has been carried out around the Scrovegni Chapel in Padua (Italy), in order to give new insights about the archaeological setting of the area. The survey made use of nearly 1500 vertical nodes deployed over two rectangular grids. 38 shot positions were fired all around the two receiver patches. The fundamental mode Rayleigh wave signal is here analysed: traveltimes are directly inferred from the signal phases, and phase velocity maps are obtained using Eikonal tomography. Also surface wave amplitudes are used, to produce autospectrum gradient maps. The joint analysis of phase velocity and autospectrum gradient allowed the identification of several buried features, among which possible remains of radial walls of the adjacent Roman amphitheater, structures belonging to a medieval convent, and the root area of an eradicated tree. Finally, depth inversion of 1D dispersion curves allowed the reconstruction of a quasi-3D shear-wave velocity model.


Introduction
This file contains a detailed description of the different steps of data processing, in order to extract Rayleigh-wave traveltimes. The full processing has been performed using original Matlab codes written by Ilaria Barone. Supplementary information also includes details about model uncertainties, and a posteriori tests on the obtainable resolution. Finally, it describes the parameters used for depth inversion of dispersion curves. This operation, as already stated in the main body of the article, is performed with Dinver, an open-source tool within the GEOPSY package 1 . The description of both data processing and dispersion-curve inversion is accompanied by figures illustrating the intermediate-step results.
Text S1. Preliminary analysis In order to perform preliminary analyses on the raw dataset, it was necessary to define a unique numbering for receiver lines inside the 3D scheme. The receiver line direction was taken as parallel to the longest side of each rectangle, which is South-East/Nort-West for receiver grid 1 and South-West/Nort-East for receiver grid 2. Numbering increases from West to East: we obtain 18 receiver lines in grid 1 (lines 1 to 18) and 20 receiver lines in grid 2 (lines 19 to 38). The analysis discussed in this section refers to all receiver lines aligned with two shot points, one at each line end, which homogeneously cover the two areas of investigation. The analyses for grid 1 and 2 were performed independently. Seismograms and f-k spectra referring to symmetrical shots inside and outside the arena are shown in Figures S1a-d and Figures S2a-d, respectively. A manual picking of each spectrum was performed to retrieve the dispersion relation for both the fundamental mode and higher modes (when clearly identifiable) for all eight shots of grid 1 and ten shots of grid 2 (Figure S1e and Figure S2e). We finally picked a dispersion curve which clearly separates fundamental modes from higher order modes for all curves (thick black line in Figure S1e and Figure S2e). This curve will serve as an LMO velocity versus frequency relation in the 3D processing (see below).

Text S2. Linear moveout correction
The linear moveout correction consists in a time shift of the recorded traces. The amount of shift has a linear dependence on the source-receiver distance (offset), being ∆t = x/V LMO , where x is the offset and V LMO is the velocity used for the correction. This simple operation flattens linear events with a velocity of V LMO , while events characterized by higher velocities acquire a negative moveout. LMO is here applied to produce a separation in the spectral domain between the fundamental mode and coherent noise (including higher modes). In particular we use the two dispersion relations, one for each receiver patch, derived in the "Preliminary analysis" section (thick black curves in Figure S1e and Figure S2e) as LMO velocity curves as a function of frequency: after the correction, fundamental mode energy at the frequency of analysis projects in the first frequency-wavenumber quadrant and higher modes in the second quadrant of the f-k plot ( Figure S3). It is important to stress how the LMO correction does not alter the position of the backscattered energy in the spectral domain, which remains in the second quadrant.

Text S3. Pseudo-2D f-k filter
Considering to the 3D nature of the acquisition and the regularity of the receiver geometry, we applied a pseudo-2D f-k filter on azimuthal sectors of 5 degrees from the source position ( Figure S4). For each slice, the true sourcereceiver distance is measured. Then, seismic traces are linearly interpolated to regularly sampled offsets and the f-k spectrum is computed. The filter sets spectral amplitudes to zero in the positive-frequency/negative-wavenumber quadrant. Finally, the filtered seismogram is retrieved through a double inverse Fourier transform and linear interpolation to the original locations. Azimuthal sectors with less than 10 traces are discarded from the analysis, as well as the near-offset region, that we consider to correspond to half wavelength. Figures S3e-f show the effect of the filter, both in the space-time and frequency-wavenumber domains: note how both higher modes and backscattering are effectively removed, while the fundamental mode energy is preserved. The effectiveness of the adopted filter is also confirmed by its ability to remove phase and amplitude oscillations, which are here mainly due to backscattering effects (see the effect of the filter on the signal phases in Figures  S5b-c). If Barone et al. 2 described this phenomenon as the result of the interference between two different modes of propagation, the periodic patterns observed here in the raw signal phases and amplitudes must be attributed to the interference between the incoming surface wave and the reflected wave from the arena middle wall. Under these conditions the oscillating period ∆X is: where k i and k r = −k i are the wavenumbers of the incoming and reflected waves. After the application of the filter,

2/18
no periodic behaviour is observed anymore, i.e. only one coherent signal (the incoming wave) and one mode of propagation are preserved.

Text S4. 2D phase unwrapping
The LMO correction already removes most of phase jumps related to the 2π periodicity of the phase. However, the LMO velocities used for the correction overestimate fundamental mode velocities, and this results in an undercorrection. Moreover, the independent application of the filter for different azimuthal sectors could introduce artifacts that translate into phase jumps at the borders between sectors. For the reasons above, residual phase unwrapping is required. We adopted a 2D unwrapping algorithm consisting of an iterative procedure. At each iteration, a phase-quality factor is computed as the sum of the variances of the phase derivatives in the X and Y directions in a close neighborhood of each point. The point with the best quality (lower variance) is selected, phase jumps larger than π are compensated in four directions (X, -X, Y, -Y) and the point is marked as "unwrapped". The unwrapping stops when all data points are unwrapped. After unwrapping, phases are de-LMO corrected and the relative traveltime maps are obtained using equation 1 of the main document. Figure S5 summarizes the different processing steps, showing their effect on surface wave phases and traveltimes.

Text S5. Phase velocity maps and autospectrum gradient maps
For the sake of completeness, we show here the phase velocity maps and the autospectrum gradient maps obtained for all frequencies of analysis ( Figures S6, S7, S8, S9).

Text S6. Phase velocity uncertainties
Eikonal tomography estimates phase velocity uncertainties as velocity standard deviations over the different shots. S10 shows the phase velocities for grid 1 and 18 Hz obtained from shot 1, the final phase velocities for the same frequency, obtained as the average of phase velocity maps for shots 1 to 19, the associated uncertainties and the number of measurements used in the analysis. It is important to stress that Figure S10b is smoother compared to Figure S10a, the latter being statistically unreliable. However, the two figures show the same major patterns, with a low velocity area in the centre. The number of measurements used for the analysis has a distribution that clearly reflects the excluded near-offset region and the side regions discarded by the f-k filter (visible also in Figure S10a), with the maximum coverage in the centre. As a consequence, standard deviations are also lower in the centre of the area. Finally, Figure S10e shows the relative error distribution, which is centered around 10% and rarely exceeds 20%, as stated in the main body of this paper ("Results" section).

Text S7. A posteriori resolution tests
The synthetic tests described in this chapter want to investigate the resolution obtainable from the application of Eikonal Tomography, given the acquisition layout used in the present study and the true (measured) uncertainties on the retrieved phase slownesses, in order to validate the obtained results. Six different synthetic models (phase velocity maps) were used: models in Figure S11a, c, e are characterized by a constant velocity of 200 m/s and an overlapping checkerboard pattern with 20% velocity variations and anomalies of size 12 m, 6 m and 3 m, respectively. Models shown in Figure S12a, c, e differ from the previous only for having lower velocity variations (10%). Both the checkerboard models and the synthetic traveltimes were generated with FMST 3 , using the same position and number of shots/receivers of the true data. We applied Eikonal Tomography to synthetic data, and perturbed single-shot slownesses with normally distributed random errors whose standard deviation is equal to the average slowness errors computed on real data (0.0008 s/m). The retrieved phase velocity maps are shown in Figure S11b, d, f and S12b, d, f. Tests show how the detectability of small scale structures is strongly influenced both by the grid spacing and by the lateral velocity contrast characterizing the anomaly, and are in a very good agreement with real data. For example, the anomaly given by the tree eradication observed in the survey area outside the amphitheatre is about 6 m wide and has a velocity contrast higher than 20%, which correspond to the model in Figure S11c. Another example is the low velocity area images at the centre of the survey area inside the amphitheatre, which is about 12 m wide and has a velocity contrast higher than 10%, corresponding to the model in Figure S12a. Despite the high level of noise, the synthetic tests demonstrate a high degree of detectability for both anomalies. These considerations on resolution based on synthetic models do not take into account possible limitations given by the frequency content of our data.

3/18
Text S8. Depth inversion We adopted a five-layer model (including a halfspace) with constant P-and S-velocities and density inside each layer. The ranges of variability of layer thicknesses and S-velocities are shown in Table S1. In particular, S-velocity ranges, depending on the average phase velocities and associated standard deviations at the different frequencies ( Figure S13), are different for grid 1 and grid 2. As for P-wave velocities, they were allowed to vary between twice the minimum V s and three times the maximum V s of each layer. Finally, the density has been kept constant for all simulations (Table S1). For each dispersion curve inversion the algorithm generated 20000 models, and computed the associated average phase velocity misfit over frequencies. The final model is the average of the best (minimum-misfit) 1000 models ( Figure S14).     1 ·V r (V r being the average phase velocity), as a function of the approximated depth, obtained as z = λ /2.5 = V r /(2.5 · f ) (λ being the average wavelength), for grid 1 and grid 2. The curves in (b) were used to define the model parameters in Table S1.